Stability analysis of a multiscale model of cell cycle dynamics coupled with quiescent and proliferating cell populations

In this paper, we perform a mathematical analysis of our proposed nonlinear, multiscale mathematical model of physiologically structured quiescent and proliferating cell populations at the macroscale and cell-cycle proteins at the microscale. Cell cycle dynamics (microscale) are driven by growth factors derived from the total cell population of quiescent and proliferating cells. Cell-cycle protein concentrations, on the other hand, determine the rates of transition between the two subpopulations. Our model demonstrates the underlying impact of cell cycle dynamics on the evolution of cell population in a tissue. We study the model’s well-posedness, derive steady-state solutions, and find sufficient conditions for the stability of steady-state solutions using semigroup and spectral theory. Finally, we performed numerical simulations to see how the parameters affect the model’s nonlinear dynamics.


Introduction
One of the cornerstones in understanding human tumor growth is mammalian cell division patterns. Many researchers have been drawn to it, and it has been the subject of extensive research for decades. Most theoretical research works explore the life cycle by utilizing agestructured frameworks. Some examples of age-structured growth models include epidemic [1][2][3], microscopic virus [4,5] and and cell population [6][7][8][9] models. However, the underlying molecular intricacies of a tissue necessitate a more comprehensive modeling framework comprising special molecular and cellular interactions.
In any living tissue, the dividing cells can be classified into quiescent and proliferating compartments. Proliferating cells divide by going through various stages in cell-cycle (G 1 , S, G 2 , M). Quiescent cells, on the other hand, do not grow or proliferate; instead, they move from the proliferative compartment to the G 0 phase and remain there until differentiation or apoptosis. For tissue homeostasis to be preserved, cells must be able to switch between the quiescent and proliferative phases. However, the transitioning between the two compartments relies on signaling molecules, which are known as growth or anti-growth factors [10]. Proliferating cells grow within a tumor cell population until the tumor is active and malignant. Besides, the total a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 [9], cell population only in proliferating phase [11,12]. Finally, cell population dynamics involving both quiescent and proliferating phases [8,[13][14][15][16][17]. Nevertheless, the influence of molecular interactions at the subcellular level on balance between proliferative and quiescent phases has not been studied. The main objective of this paper is to formulate a multiscale model by employing mathematical tools which can also encompass the heterogeneity of a complex system lying at the sub-cellular level. Therefore, we primarily focus on two predominant scales, i.e., macroscale (population dynamics) and microscale (cell-cycle dynamics), and define the coupling between these two time and age varying scales. Age refers to the time elapsed since last division, [12,18]. Note that in addition to the physical time variable (denoted by t), age-structured models introduce the age variable (denoted as a) which has rather a physiological character. The concept of "cell age" characterizes the biological variability within a proliferating cell population. Partial differential equations (PDEs) are used to simulate cell populations in the quiescent and proliferative stages at the macroscale. While ordinary differential equations are used to predict sub-cellular protein interactions related to cell-cycle dynamics (ODEs). Finally, through feedback in both directions, the two scales are connected. As mentioned earlier, proliferating cells represent a complete cycle of cell division (G 1 ,S,G 2 , M). Cells in the early proliferating phase, known as G 1 , can transition to the quiescent phase till they reach the restriction point (R). However, depending on the concentration of G 1 phase cyclin protein (x 1 ), the cells transit to S phase from late G 1 phase. It is also clear that restriction point (R) splits the cells in the G 1 -phase in two parts such that the cells become quiescent before R but can no longer avoid division once R is passed, [19,20]. In quiescent phase, cells do not divide or grow, but they continue to perform their other cellular functions. A bidirectional cell transitioning between quiescence and proliferation phases plays an essential role in tissue homeostasis, and it is regulated by extracellular environmental conditions [10]. In tumoral tissue, the balance in the bidirectional transition is disturbed, and cells may grow unconditionally [21]. Recent experiments have also revealed that cyclins are the most significant regulatory molecules for changes in cell-cycle phase, [22]. As a result, we use a crucial aspect in the dynamics of cell-cycle (i.e., from G 1 −S phase transition) to predict the evolution of a transitional balance between quiescent and proliferating subpopulations, that is essential to maintain homeostasis.
A variety of proteins are expressed at the microscale, which play an essential role in the sequential transition between different phases of cell-cycle. The complex network of protein interactions in the cell-cycle has been mathematically described using ODEs and simulated by several authors, including [23][24][25][26][27][28] and references therein. However, for simplicity, we consider only four proteins (i.e., Cyclin D À CDK4=6, p21, E2F, and Rb) from the network of proteins which participate in the cell-cycle dynamics. These proteins are chosen because they are primarily engaged in Cyclin D's activity and the progression of cells to the S from G 1 phase. The motivation stems from experimental results, which have shown that Cyclin D regulates the transition between the G 0 and G 1 phase, see [29][30][31]. Furthermore, when Cyclin D is overexpressed, cells in the proliferative phase commit to cell division, and when Cyclin D is underexpressed, cells enter a quiescent phase. It should be noted that these molecular interactions are assumed to occur in a fast growing population of cells and not in a single cell. Moreover, we assume averaged concentrations of these proteins in proliferating and quiescent cell subpopulations without considering cell to cell variability. In the sequel, we provide the biological relevance of cell-cycle proteins. The advancement in the cell-cycle is regulated by cyclin proteins (structural protein) and their cyclin-dependent kinase (CDK) inhibitors. There is a specific Cyclin/CDK complex for every phase in the cell-cycle. When micro-environment of a cell has enough growth signals, it initiates a cell-cycle that spans the activities of phase-specific complexes of cyclin protein and their catalytic partners CDK. Cyclin D activates during the G 1 phase and is induced merely by growth factors, [32]. When there are no growth factors, the concentration of Cyclin D declines, and the cell does not start the cycle. Growth-factors attach to particular receptors located on the external cytoplasmic membrane of the cell, which activates intra-cellular signaling pathways (i.e., Raf/Map/Ras kinase), which ultimately leads to the synthesis of Cyclin D (see [33][34][35], for more details). Cyclin D makes an active complex with CDK4=6 with a maximum synthesis rate. This complex can then trigger the activation of transcription factor E2F by phosphorylating its inhibitor retinoblastoma protein Rb. Resultantly, the transcription factor E2F is accumulated and activates the other essential cyclins involved in the cell-cycle.
To summarize, we develop a multiscale model to primarily address the concerns relevant to impairment in cell transitioning between quiescent and proliferating compartments, which results in unlimited tumor growth, and whether Cyclin D is responsible for the deregulation of cells transitioning between quiescent and proliferating compartments.
The layout of the paper is given in the sequel: Section 2 delves into the depth of multiscale mathematical modeling of quiescent and proliferating cell populations associated with cellcycle dynamics. In Section 3, we first demonstrate the uniqueness and existence of non-negative solutions using semigroup and spectral theory from functional analysis to confirm that the governing model equations are well-posed. Next, in Section 4, we first derive steady-state solutions and then obtain spectral criteria for local stability for steady-state solutions in a way that in the sense that if the growth bound of the linearised semigroup is negative, the steady-state solution is the locally asymptotically stable, and if growth bound is positive, the steady-state solution is unstable. Finally, Sections 5 and 6 contain the discussion about the results and final conclusion of the paper, respectively.

Age-structured model
The cell populations in quiescent and proliferating compartments are described by transport PDEs (partial differential equations) of nonlinear hyperbolic type, which characterize the density distribution of the cells concerning physiological age a and time t. In the quiescent phase, the cell density q(a, t) is given by @ @t qða; tÞ ¼ wða; x 1 Þpða; tÞ À ðgðNÞ þ m q ðaÞÞqða; tÞ; ð1Þ where the first term χ(a, x 1 )p(a, t) is the inflow from the proliferating cells at the rate χ(a, x 1 ), which is further regulated by a microscale variable, namely the age-specific concentration of Cyclin complex x 1 . The detail of the microscale variables is presented later in this section. The next term refers to the loss in quiescent cell density caused by either returning to cell division with the rate γ(N) in the proliferating phase or by cell death as a result of apoptosis (or necrosis), as depicted by death rate μ q (a). The total number cell population in both phases is represented by N(t), which is defined in Eq (3). The cells in the quiescent phase do not age (or in other words, the cells halt their age), therefore in Eq (1), the convection term concerning physiological age a is not present. In the proliferating phase, the cell number density represented by p(a, t) reads @ @t pða; tÞ þ @ @a ðgðaÞpða; tÞÞ ¼ gðNÞqða; tÞ À ðtðaÞ þ wða; where g(a) stands for the rate of evolution of a cell-cycle. The first term on the right γ(N)q(a, t) denotes the transition from the quiescent to the proliferating cells. The following term τ(a)p(a, t) symbolizes the number of cells that complete cell division at some age of the proliferating phase, whereas the cells that are moving to the quiescent phase from proliferating phase without having undergone division are given by the term χ(a, x 1 )p(a, t). Finally, the decrement in proliferating cell density because of apoptosis/necrosis is described by the death rate μ p (a). The cell population, N(t), defined as the sum of all cells in the quiescent and proliferating phases across all ages, is given as: where a ? is the maximal age of the cells. The initial conditions are defined as: The boundary condition is given as follows: for t > 0, where the number 2 appears because of the two newborn cells, which begin in the proliferating phase with age 0.The function, which defines the number of cells switching from quiescent to proliferating phase, γ(N), takes the form of monotone decreasing Hill function of N: where ν defines the maximal rate of cell transitioning from quiescent to proliferating population (e.g., when there are no cells, i.e., N = 0), κ is the Hill coefficient and θ characterises the entire cell population reaching the half maximum of ν. It means that the percentage of quiescent cells which enter the proliferative phase again declines to zero when the cell population rises, thus depicting density inhibition. The usage of the Hill function is motivated here to describe nonlinear and saturable mechanisms between the total cell population and the transition rate, see [36]. The number of cells that complete the division at some age in the proliferation phase are represented by function τ(a). The age a regulates the function τ(a) in such a way that it is almost zero until a minimum age of cells, and then it increases until the age a � : where ρ 1 is the maximum proliferation rate,ρ 2 is the age at which the half-maximum effect is achieved, and the exponent γ 1 is the Hill coefficient. Next we define the rate at which the cells leave the proliferating phase and become quiescent is given by the relation in (8). It depends on both age a and the amount of Cyclin D À CDK4=6 complex x 1 : The function χ(a, x 1 ) determines the number of cells that do not divide because of growthinhibiting factors. Age dependence in χ is motivated because the cells transit from the proliferating to quiescent phase only until a certain age that specifies a restriction point R in the cell-cycle (G 1 −S phase transition). However, until the restriction point, the concentration of Cyclin complex x 1 must be under a certain threshold to allow cells to leave the proliferating phase. In Eq (8), γ 2 and γ 3 are Hill coefficient, σ 2 and σ 3 represent the concentration of Cyclin D À CDK4=6 complex x 1 and age a, respectively, and after γ 2 and γ 3 , the rate function χ asymptotically decreases to zero and thus avoiding transition of cells to quiescent phase. It indicates that at age σ 3 , cells are inevitably devoted to entering the proliferation compartment. Lastly, σ 2 is the limit for the concentration of Cyclin' complex, which determines R, the restriction point.
In the process of cell-signaling, cell growth is regulated by the proteins called cytokine and other proliferation governing factors, [37]. Cytokines proteins attach to their special receptors, thus activating signal transduction pathways, [38]. As per different studies, it is evident that the number of cells can be kept in balance by cytokine signals, which depend on the total cell population [39]. For detailed explanation concerning the dynamics of cytokine signals, please see [40,41]. After quasi-steady-state approximation, the number of growth factors g f stemming from the total cell number N is given as, indicating maximum intensity, i.e., g f = 1, for small cell density and effectively zero intensity for large cell densities.

Cell cycle model
As previously stated, we consider only four microscale states (proteins) in the cell-cycle model, which are plausible enough to incorporate reversible transition between quiescent and proliferating phase. We utilise the kinetics of Michaelis-Menten to describe the chemical reactions with enzymes and substrates from the cell-cycle, which are briefly described in the sequel. Cyclin D protein makes a complex with its catalytic partner CDK4-6 when there are sufficient growth factors. After the formation of Cyclin D À CDK4=6 complex, it phosphorylates other proteins from the cell-cycle, which are critical to advancement in the first grwoth phase of the cell-cycle and crossing the restriction point R, [29,42]. More precisely, the Cyclin D À CDK4=6 complex phosphorylates the retinoblastoma protein Rb to inactivate it and thus release the transcription factor E2F, which in result activates many growth promoting signals to progress the cell-cycle. p21, which inhibits CDK, regulates the cell-cycle by hindering the functions of the several CDK proteins. The description of proteins is given in the Table 1. We consider the evolution of cell-cycle proteins in a single-cell whose dynamics is representative of the behavior of all cells in a population. We consider that all cells behave identical and thus one ode model with similar parameters for all cells in a population represents the microscale of underlying cell-cycle dynamics. We further postulate that our representative cell in the microscale completes division at some age a ? , while, of course, our model accounts for the cells with shorter cycles at the macroscale via function τ(a). The following ODE system describes the cell-cycle dynamics, [43]: In Eq (10a), the first term on the right-hand side describes the synthesis of Cyclin D/CDK 4-6 complex induced by the growth factors g f . The last two terms describe the binding of Cyclin D/ CDK 4-6 complex with tumor suppressor protein p21 and the degradation rate of Cyclin D/ CDK 4-6 complex, which is induced by other cell cycle proteins, for example, Cyclin E, respectively. In Eq (10b), the first term on the right-hand side describes the synthesis of transcription factors E2F induced by Cyclin D/CDK 4-6 complex. The second term denotes the decrement of E2F due to inhibition by retinoblastoma protein Rb, while the last term depicts a constant inactivation rate of E2F induced by other cell cycle proteins, for instance, Cyclin A. In the third equation (10c), the first term on the right-hand side represents the synthesis of free unphosphorylated retinoblastoma protein Rb. The second term denotes the decline in Rb by making a complex with E2F to inhibit it. The third term refers to the deactivation of Rb by phosphorylation from Cyclin D/CDK 4-6 complex and the last one to the degradation of Rb. In Eq (10d), the first and second terms represent the synthesis of p21 by ATM/ATR, TGFβ pathways and induced by E2F, respectively. The third term represents the decrement in p21 due to inhibition of Cyclin D/CDK 4-6 complex, and the last term stands for the degradation of p21. The description of the parameters involved in the cell cycle model (10a)-(10d) is described below in the Table 2. The detailed derivation of the microscale model equations is not given here; however, we suggest the interested readers to read [43] for more details. For understanding, the model simulations of above mentioned four microscale states are shown in
Proof. First, we derive the solution expression from (1) as follows: and, immediately, it follows that q(a, t) � 0 when q 0 (a) � 0 and p(a, t) is positive. Next, to derive the solution of Eq (2), we first use transformationspða; tÞ ¼ gðaÞpða; tÞ andqða; tÞ ¼ gðaÞqða; tÞ for t 2 [0, t 1 ] and a 2 [a 0 , a ? ). Then for all t 2 (0, t 1 ) and a 2 (a 0 , a ? ), we have from Eq (2): @pða; tÞ @t þ gðaÞ @pða; tÞ @a ¼ gðNðtÞÞqða; tÞ À ðtðaÞ þ wðx 1 ðaÞ; aÞ þ m p ðaÞÞpða; tÞ; ð16Þ Following that, we utilize the parameter transform given in Lemma 3.1 [40] in order to eliminate the term g(a) and introduce η as a new age variable for both p and q. We obtain @ @Zp To determine the explicit relation ofpðaðZÞ; tÞ, employ the method of characteristics (MOC). We suppose thatpðaðZÞ; tÞ is characterized by an ordinary differential equation along the curve ðaðc 1 ðyÞÞ; c 2 ðyÞÞ ¼ cðyÞ, then where c 1 ; c 2 2 R are constants. Then, it follows We can now writep using an ODE (18) so that pðað y þ c 1 Þ; y þ c 2 Þ¼pðaðc 1 ðyÞÞ; c 2 ðyÞÞ ¼ zðyÞ We may now utilize the characteristic solution to achieve the solution in {(a(η), t)|t 2 [0, t 1 ], η 2 [0, min(η � , t))}: We may now utilize the characteristic solution to achieve a solution in {(a(η), t)|t 2 [0, t 1 ], η 2 [t, η � )}: gðNðzÞÞqðaðz þ Z À tÞ; zÞdz þpðaðZ À tÞ; 0Þ � : This establishes the equation for g(a(η))p(a(η), t) for η > t. Thus, the final solution for g(a(η))p (a(η), t) can be written as: where, h(t − η) denotes the boundary conditionpðað0Þ; t À ZÞ. It can be seen that above relation for g(a)p(a, t) is positive for positive initial data and when g(a)q(a, t) � 0. Next, we check the positivity of coupled ODE model (10a)-(10d). Thereby, the set of ODEs are written as where f 1 , f 2 , f 3 and f 4 represent the vector fields of the corresponding microscale states x 1 -x 4 . Note that in Eq (19) can be seen that , which means that the concentration of x 1 increases more than it decreases for all ages. It is evident since growth factors are the only source of increase in the concentration of x 1 . Therefore, when growth factors are at the absolute minimum, x 1 is also at its lowest concentration, and hence the decrement (or degradation) cannot be more than the activation of complex x 1 . Since the solution to the system (10a)-(10b) is unique for each given initial condition (evident from (13) and (14)), thus for any given x 4 > 0, the solution will remain in the first quadrant. This guarantees the positivity of solution for x 1 . Next, we assume The solution to which takes the form x 2 ðaÞ ¼ x 0 2 e À ðk 32 x 3 ðaÞÀ k 2d Þa , which implies x 2 ðaÞ > 0 for all x 0 2 > 0 as well as for all values of x 3 ðaÞ. Thus for any given positive initial data, the solution x 2 ðaÞ is positive for all ages. In the similar fashion, we can now substitute The explicit solution cannot be computed in this case. However, the phase portrait of ðx 1 ; x 2 Þ shows that the solution trajectories point away from the axis which separate the positive and negative space for given positive initial data. In a similar way, we can also derive sufficient conditions for the positivity of the solutions for x 3 ðaÞ and x 4 ðaÞ. With this, we attain that, if z 0 2 O, z(t, z 0 ) 2 O8t > 0. Now, suppose z(t, �) = q(t, �)+ p(t, �) and death rates are identical, i.e., μ p = μ q . Then, we have from Eqs (1) and (2) where operator B generates a positive C 0 -semigroup W(t). As we know, W(t) is a nilpotent translation semigroup, it leads to z(t)(a) � q 0 (a − t) + p 0 (a − t), a > t and z(t) � 0 for t � a ? . Therefore, the mild solution z(t, z 0 ), z 0 2 O enters O 0 for t � O, and in case of z 0 2 O 0 , z(t, z 0 ) 2 O 0 , 8t � 0. Hence proved. We conclude from the above result that the norm of the local solution zðt; z 0 Þ; z 0 2 DðAÞ \ O, of (13) is defined and hence finite. As a result, we achieve the final result.
Theorem 3.3. The abstract Cauchy problem (13) has a unique global classical solution on Z with respect to the initial data z 0 2 O \ DðAÞ.
Consequently, given a positive initial data, the IVP (1) and (2) has a unique positive solution.

Existence and stability of steady-state
Here, we establish the steady-state solution of the model and sufficient conditions for the existence of the positive steady-state. First, we introduce some notations in the sequel. Let's define X as a real/complex Banach space and X ? be its dual space. The notation hF, ψi represents the value of F 2 X ? at ψ 2 X. A cone X + is defined by the following: Moreover, the dual cone, represented as X ? þ , is the subset of the dual space.

Existence of steady-states
Let � p, � q, � x 1 À � x 4 represent the steady-states of the system (1) and (2), (10a)-(10d). Then, � p, � q, � x 1 À � x 4 must satisfy these time-invariant set of ODEs: Since the ODEs of the cell-cycle model are age-dependent and with the input of growth factors at a steady-state, all cellcycle states acquire a steady-state. Therefore, to investigate the steady-states of proliferating and quiescent cell populations � pðaÞ and � qðaÞ, we do not need to solve the ODEs of the cell-cycle model explicitly. Consequently, solving the system (22) for � p and � q, we obtain � q as and after using the above relation for � q in the equation for � p yields Solving Eq (24) for � pðaÞ, yields both steady-state solutions for � pðaÞ and � qðaÞ as follows It is clear that the system defined in Eqs (1) and (2), (10a)-(10d) always admits a trivial steady-state.

Stability analysis of steady-state solutions
Next, we want to derive the stability conditions for a positive steady-state solution.Suppose qða; tÞ ¼ � q and pða; tÞ ¼ � p, 8t � 0 represent equilibrium solutions to the PDE model (1) and (2) and q � (a, t) and p � (a, t) represent the corresponding perturbation terms: qða; tÞ ¼ � q þ �q � ða; tÞ; pða; tÞ ¼ � p þ �p � ða; tÞ: Substituting the above relations in to the PDE model (1) and (2), we have where, nðtÞ≔ R a ? 0 ðp � ða; tÞ þ q � ða; tÞÞda. Then, take the derivative wrt epsilon �, leads to: Next, we formulate (25) as semilinear problem: on the Banach space X and the generator C is defined by where �ðaÞ ¼ ð� 1 ðaÞ; � 2 ðaÞÞ T 2 DðCÞ; where, D(C) is defined below: Next, the resolvent equation for operator C is considered as, Which leads to and By solving (28a), we get Which after substituting in Eq (28b) and solving gives # : Substituting ϕ 2 (a) back in Eq (29) yields where U λ and V λ are the linear operators on Banach space, given as Similarly, we rewrite ϕ 2 (a) as Let L ¼ fl 2 C j À m q ð�Þ À gð � N Þ 2 sðU l Þg, then we can say that if l 2 C n L, operators U λ and V λ are compact operators from X to L 1 (0, a ? ). This implies ϕ 1 (a) is represented by a compact operator. In a similar fashion, ϕ 2 (a) is also represented by a compact operator. Resultantly, we get that operator C has a compact resolvent which further implies that σ(C) comprises entirely of isolated eigenvalues, i.e., σ(C) = σ P (C) (see p. 187, Theorem 6.29 in [46]). From latter, we know that C n L � rðCÞ, where ρ(C) is the resolvent of C. This implies σ P (C) = σ(C) � Λ. Since U λ is a compact operator, then it leads to σ(U λ )\{0} = σ P (U λ )\{0}. Now if λ 2 Λ, there exists an eigenfunction ψ λ such that U λ ψ λ = ψ λ . Then, it is trivial to see that (ϕ 1 (a), ϕ 2 (a)) T provides an eigenvector of C for an eigenvalue λ. Then Λ � σ P (C), and finally, we can say that (30) satisfies.

Lemma 4.2. Let T(t) be the C 0 -semigroup generated by the operator C, t � 0. Then, T(t) is eventually norm continuous (ENC) and
where ω 0 (C) represents the growth bound of semigroup T(t) and s(C) denotes the spectral bound of the operator C.
Proof. First, we write the bounded operator C as: for ϕ 2 X. To prove the compactness of C, we show that for any bounded sequence ð� n Þ n2N in X, the sequence ðC� n Þ n2N has a uniformly convergent subsequence. For this we use the Arzelà-Ascoli Theorem. Thereby, we need to check that ðC� n Þ n2N is uniformly bounded and uniformly equicontinuous. For the boundedness, note that since we assumed that ð� n Þ n2N is bounded, we have kC� n k 1 � kCk k� n k 1 � kCk sup n2N k� n k 1 ; proving that ðC� n Þ n2N is also bounded. Next, for the uniform equicontinuity, consider Z R jðC�Þða þ hÞ À ðC�ÞðaÞjda ¼ where kðaÞ ¼ À @ @a gðaÞ À tðaÞ À wða; � x 1 Þ À m p ðaÞ. It follows that ðC� n Þ n2N is equicontinuous. Thus, by the Arzelà-Ascoli Theorem, the sequence ðC� n Þ n2N has a uniformly convergent subsequence, and therefore, C is compact which implies T is ENC semigroup. As we know that the spectral mapping theorem applies to ENC semigroup, we get the spectral determined growth condition, i.e., ω 0 (C) = s(C), thus we obtain (33).
Next, to study the stability of equilibrium states, we need to find that the dominant singular point, i.e., the element of set Λ which has the largest real part. Then utilizing (30) and (33), we can find the growth bound of semigroup T. Lemma 4.3. The operator U λ , l 2 R is nonsupporting with respect to X + and lim l!þ1 holds. Proof. It can be seen from (31) and (32) that the operator U l ; l 2 R is strictly positive. Now, in order to show non-supporting property of U l ; l 2 R, we can easily verify the inequality where the linear function f λ , is given as Thereby, it leads us to U nþ1 l c � hf l ; cihf l ; ci n c; 8n: Since f λ is strictly positive and the constant function c = 1 is a quasi-interior point of L 1 (0, a ? ), it leads to hF; U n l i > 0 for every pair ψ 2 X + \{0}, F 2 X � þ n f0g. Then U l ; l 2 R is nonsupporting. Following that, we utilise (35) and take duality pairing with the eigenfunctional F λ of U λ which corresponds to r(U λ ), then we get rðU l ÞhF l ; ci � hF l ; eihf l ; ci: Suppose ψ = c, we obtain an inequality r (U λ ) � hf λ , ci, where It follows that By using the positivity of gð � N Þ; m p ; m q ; w and τ, we conclude the following lim l!þ1 rðU l Þ ¼ 0: Hence proved. Preceding Lemma concludes that λ ! r(U λ ) is a decreasing function of l 2 R. Furthermore, if l 2 R so that r(U λ ) = 1, then λ 2 Λ since r(U λ ) 2 σ P (U λ ). From the monotonicity of r (U λ ) and (34), the following holds.

Results and discussion
This section presents some of the model simulation to understand the evolution of both subpopulations in relation with the cell-cycle dynamics. Table 3 shows the model parameters employed in the simulations. Most of the parameters are used from the literature, however the rest of them were unknown and we used some arbitrary values which are selected either likely to be biologically relevant or by using a range of values so that our numerical simulations exhibit the expected behavior of the model under implied assumptions. Please note that, we do not consider a specific organ while choosing parameters, instead we consider the parameters for stem cell lines and for experimental validation, one must identify these parameters for other cell lines. Furthermore, we assume that the maximal age a ? of the cells is 50. The spatial step size and the time step is used as Δa = 0.5 and Δt = 0.02, respectively. Additionally, for the sake of clarity, we also assume unit speed, i.e., g(a) = 1. Three case studies will be discussed in the following section.

Local stability of a non-trivial steady-state solution
Here, we investigate local stability of the non-trivial equilibrium solution. The parameter values used are μ p = μ q = 0.0014. We take γ(N) = 6.8964 × 10 −6 and ρ 1 = 1.0. The initial conditions are assumed as qða; 0Þ ¼ pða; 0Þ ¼ k 0 ffi ffi ffi ffi ffi ffiffi   (9) which depends on total cell population is an ideal representation of growth factors which entirely relies on total number of cells and ignores various other possible scenarios, for instance, availability of nutrients, PH level, oxygen concentration etc. Furthermore, the gamma function γ(N) which determines the cell transitions to proliferating from quiescent cells, is shown in Fig 4(c). It represents an inverse relation to total cell population and declines to a very low number when the respective cell populations attain a steady-state. In terms of feedback from cell-cycle to population level, merely Cyclin D À CDK4=6 complex' concentration is taken into account. It mainly influences the transition rate χ(a, x 1 ) from proliferating to quiescent cells. It is evident from the distribution of proliferating cells in Fig 3 that new cells are entering proliferating phase at age a = 0 and after 20 hours of aging, cells start leaving the proliferating phase depending on their cycle length and concentration of the complex Cyclin D À CDK4=6. However, quiescent cells q(a, t) are accumulating in the early proliferating phase which do not achieve certain threshold of Cyclin D À CDK4=6 concentration to pass through restriction point in the cell-cycle.

Local stability of the trivial solution
Next, we investigate the local stability of the trivial equilibrium solution. Thereby, we choose the death rates to be constants and μ p = μ q = 0.020. Moreover, we take ρ 1 = 0.20 and ν = 0.1.
We used the initial conditions as follows pða; 0Þ ¼ qða; 0Þ ¼ k 0 ffi ffi ffi ffi ffi ffiffi  Fig 5(a). The trivial steady-state is achieved until 2500 hours and cell population declines to zero. The growth factors, on the other hand, reach to their maximum value 1 and retain that value throughout due to very low cell number. The gamma function also attains its maximum with time.

Instability in the solution
Eventually, in Fig 6, we explore the parameters which lead to an instability in the solution.
The proposed model, in general, displays very robust dynamics because of the closed feedback-loops. Nonetheless, transitioning function like χ(a, x 1 ) display some sensitivity to the fluctuations in cell-cycle states. In order to pose a scenario to see abnormal cell-cycle behavior, we altered a parameter k gf = 0.0001, which reflects a situation in which synthesis of Cyclin D À CDK4=6 complex is somewhat altered due to the change in the influence of growth factors on it. More precisely, by altering the parameter k gf , we induce delays in the oscillation of Cyclin D À CDK4=6 complex. Consequently, the cell number increases exponentially. Rest of the parameters are similar to other case studies explained before. Total cell population is plotted in Fig 6(a) which increases rapidly when there are abundant growth factors, see Fig 6(b). Finally, the transition function γ(N) is declining with time, as expected.  The proposed model does have some limitations also. The model, for example, excludes cell-to-cell variability, which is an important aspect to capture noise and heterogeneity from the cellular level. The feedback model, which includes growth factors, is relatively simple, and activation of the Cyclin D complex can only be characterized by taking into account all signaling pathways. Furthermore, at the microscale, the cell-cycle model is confined only to Cyclin D and the proteins in direct interaction with it; nevertheless, multiple additional proteins can control this network in various situations. Finally, while the Cyclin D complex and its inhibitor CDK4=6 plays a crucial part in the G 1 to S phase transition, the other restriction point in the S phase for detecting DNA damage has been overlooked.

Conclusion
This study presents non-linear, multiscale modeling of physiologically-structured quiescent and proliferating cells in relation to cell-cycle dynamics, which play an essential part in committing a cell to irreversible cell-division process. We assume reversible transitioning from quiescent to proliferating cells and vice versa and, additionally, a feedback in both directions which maintains the homeostasis. We checked the wellposedness of the model, derive non-trivial equilibrium solutions and find spectral criteria for local stability in the sense that if the growth bound of the linearised semigroup is negative, the steadystate solution is the locally asymptotically stable, and if growth bound is positive, the steady-state solution is unstable. We also performed numerical simulations to study the behavior of the proposed model, and, in this regard, we studied three scenarios with some variation in the parameters. The first scenario explains the steady-state behavior of the model in any healthy person under normal conditions. The second scenario relates to a trivial steady-state where, hypothetically, the decline in cell number density is more than the rise due to newborn cells. Finally, in the third case study, we investigate the impact of Cyclin D À CDK4=6 complex on the transition between two sub-populations. It turns out that any fluctuations in synthesis and degradation of Cyclin D À CDK4=6 complex can result in an abnormal growth in cell number, thus leading to cancer. Moreover, it shows that the Cyclin complex plays a vital role in the reversible transition between the two subpopulations.
For possible future extensions of this work, we intend to extend the current modeling framework in an optimal control problem setting and to perform a thorough sensitivity analysis of the involved parameters.
In the similar fashion, F 1 and F 2 can be shown Fréchet differentiable in both X and Y.